### Making Sense of Human Rights Diplomacy: Evidence from a US Campaign To Free Political Prisoners
## Myrick & Weinstein, International Organization

### Replication File

# Load packages

library(foreign)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(lubridate)
library(stringr)
library(ebal)
library(Matching)
library(MASS)
library(rgenoud)
library(survminer)
library(survival)
library(lattice)
library(stargazer)
library(openxlsx)
library(tidyverse)
library(ggthemes)
library(grid)
library(gridExtra)
library(msm)
library(estimatr)
library(lmtest)
library(coefplot)
library(sandwich)
library(multiwayvcov)
library(panelView)
library(lfe)
library(dfadjust)
library(openxlsx)
library(questionr)

# Set working directory

setwd("~/Dropbox/ATI Free the 20/Replication/")

# Figure 2: Short List vs. Long List Outcomes -----------------------------

# Create a survival curve for F20 post-Sept. 2015
s = read.csv("f20_prisoners_data.csv", stringsAsFactors = FALSE)
s$start = s$t1
s$start = parse_date_time(s$start, "m!/d!/y!")

s$end = s$t2
s$end = parse_date_time(s$end, "m!/d!/y!")

s$launch = "9/01/15"
s$launch = parse_date_time(s$launch, "m!/d!/y!")
s$months_since_arrest = round(difftime(s$end, s$start, units = "days")/30)
s$months_since_f20 = round(difftime(s$end, s$launch, units = "days")/30)

s$`Short List` = s$f20

# T-Test for Release Outcome (Difference in Means)

short = subset(s, f20==1)
long = subset(s, f20==0)

table(short$released) # 17 of 19 released on "short list" (89%)
table(long$released) # 11 of 21 released on "long list" (52%)

t.test(short$released, long$released) 
t.test(short$released, long$released[long$released_before_launch==0])

# Figure 2 (K-M Plot of Release Outcomes for "Short" vs. "Long" List)

s2 = subset(s, released_before_launch == 0)

sfit2 = survfit(Surv(s2$months_since_f20, s2$released)~s2$f20, data=s2)
survplot2 = ggsurvplot(sfit2, data = s2, risk.table = FALSE, conf.int = TRUE, 
                       palette=c("gray", "black"),
                       ggtheme = theme_minimal(base_size=30), 
                       pval = TRUE, tables.theme = clean_theme(),
                       legend.title = "", legend.labs = c("Long List", "Short List"))
survplot2 + ylab("Probability") + xlab("Time Since #Freethe20 (Months)") + 
    ggtitle("") 

# Figure 3: Count of Females by Country ---------------------------------------------

# Prep UA Data

f = read.csv("amnesty_ua_data.csv", stringsAsFactors = FALSE)
f = subset(f, select = c(ua_no, id_no, date, female, name, title, country, 13:54))
f = f[-which(is.na(f$female)),] #remove unnamed individuals

# Plot Females vs. Males in UA Data

g1 = ggplot(f, aes(as.factor(country)))
g1 + theme_bw() + geom_bar(aes(fill=as.factor(female))) + 
    theme(axis.text.x = element_text(angle=65, vjust=0.5, size = 20), 
          axis.text.y = element_text(size = 15),
          axis.title = element_text(size = 25),
          plot.title = element_text(size = 25),
          legend.text = element_text(size = 20)) + 
    labs(title="Women in Amnesty International Urgent Actions (2000-2015)", y = "Count", x = "Country") + 
    scale_fill_manual(name="", labels=c("Men", "Women"), values = c("gray", "black"))


# Code Control Variables --------------------------------------------------

# Any Medical Concerns

f$any_medical = ifelse(f$access_to_medical == 1 | f$denial_medical == 1 |
                           f$health_concern == 1 | f$medical_concern == 1 |
                           f$no_medical == 1, 1, 0)

# Any Torture Concerns

f$any_torture = ifelse(f$fear_of_torture == 1 | f$torture == 1, 1, 0)

# Any Death Penalty Concerns

f$any_death_penalty = ifelse(f$death_penalty == 1 | f$fear_execution == 1 | 
                                 f$imminent_execution == 1, 1, 0)

# Reported as a Prisoner of Conscience

f$any_poc = ifelse(f$poc== 1 | f$poss_poc == 1, 1, 0)

# Any Trial or Legal Concerns

f$any_trial_legal = ifelse(f$fear_unfair_trial == 1 | f$legal_concern == 1 |
                               f$unfair_trial == 1 | f$detention_notrial == 1 |
                               f$poss_unfair_trial == 1, 1, 0)

# Count number of appearances in dataset by name

f = f %>% add_count(name)
f$count_by_name = f$n

# Recode date and year variables

f$date = parse_date_time(f$date, "d!-Om!*-y!*")
f$year = year(f$date)

# Code earliest and latest appearance in dataset by name

f = f %>% group_by(name) %>% 
    mutate(earliest_date = min(date, na.rm = TRUE)) %>% ungroup()

f = f %>% group_by(name) %>% 
    mutate(latest_date = max(date, na.rm = TRUE)) %>% ungroup()

View(subset(f, select = c(name, date, earliest_date, latest_date)))

# Reshape Data by Name

fn = subset(f, select = c(female, name,
                          country, free20, any_poc, 
                          any_torture, any_death_penalty, 
                          any_trial_legal, any_medical, 
                          count_by_name, earliest_date, latest_date))

fn = fn %>% group_by(name) %>% mutate(any_poc = max(any_poc, na.rm = TRUE),
                                      any_torture = max(any_torture, na.rm = TRUE),
                                      any_death_penalty = max(any_death_penalty, na.rm = TRUE),
                                      any_trial_legal = max(any_trial_legal, na.rm = TRUE),
                                      any_medical = max(any_medical, na.rm = TRUE)) %>%
    ungroup()

fn = fn[order(fn$country, fn$earliest_date),]
fn = fn[!duplicated(fn), ]

fn$year_entered = as.numeric(str_sub(as.character(fn$earliest_date), 1, 4))

# Create Country Dummies

fn$azerbaijan = ifelse(fn$country == "Azerbajian", 1, 0)
fn$china = ifelse(fn$country == "China", 1, 0)
fn$egypt = ifelse(fn$country == "Egypt", 1, 0)
fn$eritrea = ifelse(fn$country == "Eritrea", 1, 0)
fn$ethiopia = ifelse(fn$country == "Ethiopia", 1, 0)
fn$iran = ifelse(fn$country == "Iran", 1, 0)
fn$myanmar = ifelse(fn$country == "Myanmar", 1, 0)
fn$russia = ifelse(fn$country == "Russia", 1, 0)
fn$syria = ifelse(fn$country == "Syria", 1, 0)
fn$uzbekistan = ifelse(fn$country == "Uzbekistan", 1, 0)
fn$venezuela = ifelse(fn$country == "Venezuela", 1, 0)
fn$vietnam = ifelse(fn$country == "Vietnam", 1, 0)

# Subset only to female political prisoners

fnf = subset(fn, female==1)
fnf = fnf[!duplicated(fnf), ] 

#Note: This is the base dataset for female_prisoners_data.xlsx 
#It is then used to manually code arrest/releases and check for transliteration/duplicates in the base data

#write.csv(fnf, "female_prisoners_data_base.csv")

# female_prisoners_data_base.csv contains no arrest/release information
# female_prisoners_data.xlsx contains arrest/release information and other manually coded information

# Table 3: Logistic Regressions -----------------------------------------------

c = read.xlsx("female_prisoners_data.xlsx")
c = subset(c, prisoner==1 & female==1) #remove women who are not political prisoners (manually coded, see appendix for coding rules)

# Create start and end time variables
c$start = ifelse(!is.na(c$start), paste( "01", as.character(c$start), sep = "/"), c$start)
c$start = parse_date_time(c$start, "d!/m!/y!")

c$end = ifelse(!is.na(c$end), paste( "01", as.character(c$end), sep = "/"), c$end)
c$end = ifelse(is.na(c$end) & c$not_released==1, "01/09/2018", c$end)
c$end = parse_date_time(c$end, "d!/m!/y!")

# Create censored variable (1 = censored (not released); 2 = released)
c$status = ifelse(c$not_released==1, 1, NA)
c$status = ifelse(c$not_released==0, 2, c$status)

#Create binary indicator for released
c$released = ifelse(c$status==1, 0, NA)
c$released = ifelse(c$status==2, 1, c$released)

#Create time elapsed variable
c$time = round(difftime(c$end, c$start, units = "days")/30)

#Remove cases missing start date or end date 
c = subset(c, !is.na(time))

#Create variable for F20
c$launch = "9/01/15"
c$launch = parse_date_time(c$launch, "m!/d!/y!")
c$time_since_f20 = round(difftime(c$end, c$launch, units = "days")/30)

c$cntrl1 = ifelse(c$time_since_f20 > 0, 1, 0)

#Create censoring variable if conducting analysis in Stata
c$event = ifelse(c$status==2, 1, NA)
c$event = ifelse(c$status==1, 0, c$event)

#Create regions variable
c$eurasia = ifelse(c$azerbaijan == 1 | c$uzbekistan == 1 | c$russia==1, 1, 0)
c$asia = ifelse(c$china == 1 | c$myanmar == 1 | c$vietnam==1, 1, 0)
c$mideast = ifelse(c$egypt == 1 | c$iran == 1 | c$syria==1, 1, 0)
c$americas = ifelse(c$venezuela == 1, 1, 0)
c$africa = ifelse(c$eritrea==1 | c$ethiopia == 1, 1, 0)

c$`Free the 20` = c$f20

# Logistic Regression Models: Outcome is Release=1

# Sample is all female political prisoners (Models 1-4)
# Or all female political prisoners in jail at F20 launch (Model 5)

# Table 3: Were #Freethe20 women more likely to get released? 

b1r = glm(released ~ free20, data = c)
b2r = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical, data = c)
b3r = glm(released ~ free20 + africa + asia + mideast, data = c)
b4r = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + africa + asia + mideast, data = c)
b5r = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + africa + asia + mideast, data = subset(c, time_since_f20 > 0))

stargazer(b1r, b2r, b3r, b4r, b5r, title = "Were Freethe20 women more likely to get released?", 
          covariate.labels = c("Freethe20", "Prisoner of Conscience",
                               "Any Torture Concerns", "Any Death Penalty Concerns",
                               "Any Legal Concerns", "Any Medical Concerns",
                               "Africa", "Asia", "Middle East"),
          label = "tab:logistic",
          no.space = TRUE,
          omit.stat = c("chi2"))

predict(b1r, data.frame(free20 = 1), type="response") #89
predict(b1r, data.frame(free20 = 0), type="response") #43

# Table 4: Weibull Models ----------------------------------------------------------

# Survival analysis with Weibull models
# Sample is all political prisoners (Models 1-4) or political prisoners in jail at F20 launch (Model 5)
# Outcome is time since arrest (Models 1-4) or time since F20 (Model 5)

# Table 4: Were #Freethe20 women released at a faster rate? 

w1r = survreg(Surv(time, status)~free20, data = c, dist = 'weibull')

w2r = survreg(Surv(time, status)~free20+ any_poc + any_torture + any_death_penalty + 
                 any_trial_legal + any_medical, data = c, dist = 'weibull')

w3r = survreg(Surv(time, status)~free20 + africa + asia + mideast, data = c, dist = 'weibull')

w4r = survreg(Surv(time, status)~free20 + any_poc + any_torture + any_death_penalty + 
                 any_trial_legal + any_medical + africa + asia + mideast, data = c, dist = 'weibull')

w5r = survreg(Surv(time_since_f20, status)~free20 + any_poc + any_torture + any_death_penalty + 
                 any_trial_legal + any_medical+ africa + asia + mideast, data = subset(c, time_since_f20>0), dist = 'weibull')

stargazer(w1r, w2r, w3r, w4r, w5r, title = "Survival Analysis for Female Political Prisoners", 
          covariate.labels = c("Freethe20", "Prisoner of Conscience",
                               "Any Torture Concerns", "Any Death Penalty Concerns",
                               "Any Legal Concerns", "Any Medical Concerns",
                               "Africa", "Asia", "Middle East"),
          label = "tab:weibull",
          no.space = TRUE,
          omit.stat = c("chi2"))

# Figure 4: Media Coverage ------------------------------------------------

m = read.xlsx("f20_media_data.xlsx")

m$f20_month_lag1 = ifelse(m$month==10 & m$year==2015, 1, 0)
m$f20_month_lag2 = ifelse(m$month==11 & m$year==2015, 1, 0)
m$f20_month_lead1 = ifelse(m$month==8 & m$year==2015, 1, 0)
m$f20_month_lead2 = ifelse(m$month==7 & m$year==2015, 1, 0)

# Media Coverage (Lexis Nexis Data)

# The dependent variable is equal to ln(1 + total news articles per month)  
# SEs are 1 and 2 robust standard errors (results robust to alternative SEs)

lexnex = lm_robust(media_ln ~ f20_month_lead2 + f20_month_lead1 + f20_month + 
                       f20_month_lag1 + f20_month_lag2 + arrest + 
                       release + as.factor(year) + as.factor(name_id), data = m)


lnplot = coefplot(lexnex, coefficients = c("f20_month_lead2", "f20_month_lead1", 
                                           "f20_month", "f20_month_lag1",
                                           "f20_month_lag2"),
                  innerCI=1, outerCI=2, lwdInner=2, lwdOuter=1,
                  color = "black")

lnplot + coord_flip() + theme_minimal(base_size=20) + 
    ggtitle("") + xlab("Total News Articles (ln)") + ylab("") +
    scale_y_discrete(labels=c("July 2015", "Aug 2015", 
                              "Sept 2015 \n (Free the 20)", 
                              "Oct 2015", "Nov 2015"))

# Online Search Interest (Google Trends Data)

# The dependent variable is equal to ln(1 + maximum Google Trend Value)
# SEs are 1 and 2 robust standard errors (results robust to alternative SEs)

google = lm_robust(max_goog_ln ~ f20_month_lead2 + f20_month_lead1 + f20_month + 
                       f20_month_lag1 + f20_month_lag2 + arrest + 
                       release + as.factor(year) + as.factor(name_id), data = m)

gplot = coefplot(google, coefficients = c("f20_month_lead2", "f20_month_lead1", 
                                          "f20_month", "f20_month_lag1",
                                          "f20_month_lag2"),
                 innerCI=1, outerCI=2, lwdInner=2, lwdOuter=1,
                 color = "black")

gplot + coord_flip() + theme_minimal(base_size=20) + 
    ggtitle("") + xlab("Maximum Google Trend Value (ln)") + ylab("") +
    scale_y_discrete(labels=c("July 2015", "Aug 2015", 
                              "Sept 2015 \n (Free the 20)", 
                              "Oct 2015", "Nov 2015"))

# Appendix V --------------------------------------------------------------

# Comparison of F20 Women to Other Female Political Prisoners Featured by Amnesty International

# Read in UA Arrest and Release Data

X = as.matrix(subset(c, select = c("count_by_name", "any_poc", "any_torture", 
                                   "any_death_penalty", "any_trial_legal", "any_medical")))

# Compare characteristics of F20 vs. AI women in same countries

mb = MatchBalance(free20 ~ X, data = c)

out = baltest.collect(matchbal.out=mb,var.names=colnames(X), after=FALSE)
out = subset(out, select = c(mean.Tr, mean.Co, `T pval`))

library(xtable)
xtable(round(out, 3))

# Appendix VI --------------------------------------------------------------

# Robustness to Dropping Torture Cases

b1t = glm(released ~ free20, data = subset(c, any_torture==0))
b2t = glm(released ~ free20 + any_poc + any_death_penalty + any_trial_legal + any_medical, 
          data = subset(c, any_torture==0))
b3t = glm(released ~ free20 + africa + asia + mideast, data = subset(c, any_torture==0))
b4t = glm(released ~ free20 + any_poc + any_death_penalty +
              any_trial_legal + any_medical + africa + asia + mideast, data = subset(c, any_torture==0))
b5t = glm(released ~ free20 + any_poc + any_death_penalty +
              any_trial_legal + any_medical + africa + asia + mideast, data = subset(c, any_torture==0 & time_since_f20 > 0))

stargazer(b1t, b2t, b3t, b4t, b5t, title = "Were Freethe20 women more likely to get released? (No Torture Cases)", 
          covariate.labels = c("Freethe20", "Prisoner of Conscience",
                               "Any Death Penalty Concerns",
                               "Any Legal Concerns", "Any Medical Concerns",
                               "Africa", "Asia", "Middle East"),
          label = "tab:logistic",
          no.space = TRUE,
          omit.stat = c("chi2"))

# Robustness to Dropping Death Penalty Cases

b1d = glm(released ~ free20, data = subset(c, any_death_penalty==0))
b2d = glm(released ~ free20 + any_poc + any_torture + any_trial_legal + any_medical, 
          data = subset(c, any_death_penalty==0))
b3d = glm(released ~ free20 + africa + asia + mideast, data = subset(c, any_death_penalty==0))
b4d = glm(released ~ free20 + any_poc + any_torture + 
              any_trial_legal + any_medical + africa + asia + mideast, data = subset(c, any_death_penalty==0))
b5d = glm(released ~ free20 + any_poc + any_torture + 
              any_trial_legal + any_medical + africa + asia + mideast, data = subset(c, any_death_penalty==0 & time_since_f20 > 0)) 

stargazer(b1d, b2d, b3d, b4d, b5d, title = "Were Freethe20 women more likely to get released? (No Death Penalty Cases)", 
          covariate.labels = c("Freethe20", "Prisoner of Conscience",
                               "Any Torture Concerns", 
                               "Any Legal Concerns", "Any Medical Concerns",
                               "Africa", "Asia", "Middle East"),
          label = "tab:logistic",
          no.space = TRUE,
          omit.stat = c("chi2"))

# Regression with Country Fixed Effects

b1 = glm(released ~ free20, data = c)
b2 = glm(released ~ free20 + as.factor(country), data = c) 
b3 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
           any_trial_legal + any_medical, data = c)
b4 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
           any_trial_legal + any_medical + as.factor(country), data = c)
b5 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
           any_trial_legal + any_medical, data = subset(c, time_since_f20 > 0))
b6 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
           any_trial_legal + any_medical + as.factor(country), data = subset(c, time_since_f20 > 0))

stargazer(b1, b2, b3, b4, b5, b6, title = "Were Freethe20 women more likely to get released?", 
          label = "tab:logistic",
          omit = c("country"), 
          covariate.labels = c("Freethe20", "Prisoner of Conscience",
                               "Any Torture Concerns", "Any Death Penalty Concerns",
                               "Any Legal Concerns", "Any Medical Concerns"),
          no.space = TRUE,
          omit.stat = c("chi2"),
          add.lines = list(c("Country FEs", "No", "Yes", "No", "Yes", "No", "Yes")))


# Appendix VII ------------------------------------------------------------


# Note: Sensitivity analysis in Appendix VII was conducted in Shiny app from Cinelli & Hazlett (2020)
# Source: https://carloscinelli.shinyapps.io/robustness_value/

# To compute bounding strength, use full model from Table 3 (Model 4) in manuscript

    # Covariate used for bounding: any_death penalty
    # t-value of covariate in outcome regression: -2.031
    # t-value of covariate in treatment regression: -3.166

summary(glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
        any_trial_legal + any_medical + africa + asia + mideast, data = c))

summary(glm(free20 ~ any_poc + any_torture + any_death_penalty +
                any_trial_legal + any_medical + africa + asia + mideast, data = c))

# Appendix VIII --------------------------------------------------------------

## Robustness checks for survival models

# COX PROPORTIONAL HAZARD MODELS

cox1 = coxph(Surv(time, status)~free20, data=c)

cox2 = coxph(Surv(time, status)~free20 + any_poc + any_torture + any_death_penalty + 
                 any_trial_legal + any_medical, data = c)

cox3 = coxph(Surv(time, status)~free20 + africa + asia + mideast, 
             data = c)

cox4 = coxph(Surv(time, status)~free20 + any_poc + any_torture + any_death_penalty + 
                 any_trial_legal + any_medical + africa + asia + mideast, 
             data = c)

cox4b = coxph(Surv(time_since_f20, status)~free20 + any_poc + any_torture + any_death_penalty + 
                  any_trial_legal + any_medical + africa + asia + mideast, 
              data = subset(c, cntrl1==1))

# Note: Proportional hazard assumption does not hold, which is why these models
# are not used in the main analyses

cox.zph(cox1)
cox.zph(cox2)
cox.zph(cox3)
cox.zph(cox4)
cox.zph(cox4b)

stargazer(cox1, cox2, cox3, cox4, cox4b,
          title = "Were Freethe20 women released at a faster rate? (Cox Models)", 
          covariate.labels = c("Freethe20", "Prisoner of Conscience",
                     "Any Torture Concerns", "Any Death Penalty Concerns",
                     "Any Legal Concerns", "Any Medical Concerns",
                     "Africa", "Asia", "Middle East"),
          label = "tab:cox",
          no.space = TRUE)


# LOG-NORMAL MODELS

ln1 = survreg(Surv(time, status)~free20, data = c, dist = 'lognormal')

ln2 = survreg(Surv(time, status)~free20+ any_poc + any_torture + any_death_penalty + 
                 any_trial_legal + any_medical, data = c, dist = 'lognormal')

ln3 = survreg(Surv(time, status)~free20 + africa + asia + mideast, data = c, dist = 'lognormal')

ln4 = survreg(Surv(time, status)~free20 + any_poc + any_torture + any_death_penalty + 
                 any_trial_legal + any_medical + africa + asia + mideast, data = c, dist = 'lognormal')

ln5 = survreg(Surv(time_since_f20, status)~free20 + any_poc + any_torture + any_death_penalty + 
                 any_trial_legal + any_medical + africa + asia + mideast, data = subset(c, time_since_f20>0), dist = 'lognormal')

stargazer(ln1, ln2, ln3, ln4, ln5, title = "Were Freethe20 women released at a faster rate? (Log-Normal Models)", 
          covariate.labels = c("Freethe20", "Prisoner of Conscience",
                               "Any Torture Concerns", "Any Death Penalty Concerns",
                               "Any Legal Concerns", "Any Medical Concerns",
                               "Africa", "Asia", "Middle East"),
          label = "tab:lognormal",
          no.space = TRUE)

# LOG LOGISTIC MODELS

ll1 = survreg(Surv(time, status)~free20, data = c, dist = 'loglogistic')

ll2 = survreg(Surv(time, status)~free20+ any_poc + any_torture + any_death_penalty + 
                  any_trial_legal + any_medical, data = c, dist = 'loglogistic')

ll3 = survreg(Surv(time, status)~free20 + africa + asia + mideast, data = c, dist = 'loglogistic')

ll4 = survreg(Surv(time, status)~free20 + any_poc + any_torture + any_death_penalty + 
                  any_trial_legal + any_medical + africa + asia + mideast, data = c, dist = 'loglogistic')

ll5 = survreg(Surv(time_since_f20, status)~free20 + any_poc + any_torture + any_death_penalty + 
                  any_trial_legal + any_medical + africa + asia + mideast, data = subset(c, time_since_f20>0), dist = 'loglogistic')

stargazer(ll1, ll2, ll3, ll4, ll5, title = "Were Freethe20 women released at a faster rate? (Log-Logistic Models)", 
          covariate.labels = c("Freethe20", "Prisoner of Conscience",
                               "Any Torture Concerns", "Any Death Penalty Concerns",
                               "Any Legal Concerns", "Any Medical Concerns",
                               "Africa", "Asia", "Middle East"),
          label = "tab:loglogistic",
          no.space = TRUE)

# Appendix IX -------------------------------------------------------------

library(estimatr)

## Robustness checks for media coverage models

###### Online Search Interest (Google Trend Values)

# Change DV to Google Trend Value 

g1 = lm(max_goog ~ f20_month_lead2 + f20_month_lead1 + f20_month + 
             f20_month_lag1 + f20_month_lag2 + arrest + 
             release + as.factor(year) + as.factor(name_id), data = m)

# Change DV to ln(1 + max GTV per month)  

g2 = lm(max_goog_ln ~ f20_month_lead2 + f20_month_lead1 + f20_month + 
            f20_month_lag1 + f20_month_lag2 + arrest + 
            release + as.factor(year) + as.factor(name_id), data = m)

stargazer(g1, g2, se = starprep(g1, g2),
          covariate.labels = c("July 2015", "August 2015", "September 2015 (F20)", 
                               "October 2015", "November 2015", 
                               "Arrest", "Release"),
          title = "Did Freethe20 increase online searches of the featured women?",
          label = "tab:goog", 
          omit = c("year", "name_id", "Constant"),
          no.space = TRUE)

###### Media Coverage (Lexis Nexis Data)

# Change DV to total news articles per month

lx1 = lm(media ~ f20_month_lead2 + f20_month_lead1 + f20_month + 
             f20_month_lag1 + f20_month_lag2 + arrest + 
             release + as.factor(year) + as.factor(name_id), data = m)
# Change DV to ln(1 + total news articles per month)  

lx2 = lm(media_ln ~ f20_month_lead2 + f20_month_lead1 + f20_month + 
                    f20_month_lag1 + f20_month_lag2 + arrest + 
                    release + as.factor(year) + as.factor(name_id), data = m)

stargazer(lx1, lx2, se = starprep(lx1, lx2),
          covariate.labels = c("July 2015", "August 2015", "September 2015 (F20)", 
                               "October 2015", "November 2015", 
                               "Arrest", "Release"),
          title = "Did Freethe20 increase media coverage of the featured women?",
          label = "tab:lexnex",
          omit = c("year", "name_id", "Constant"),
          no.space = TRUE)

###### Count Models: Google Trends & Lexis Nexis Data

nb1 = glm.nb(max_goog ~ f20_month_lead2 + f20_month_lead1 + f20_month + 
                   f20_month_lag1 + f20_month_lag2 + arrest + 
                   release + as.factor(year) + as.factor(name_id), data = m)

nb2 = glm.nb(media ~ f20_month_lead2 + f20_month_lead1 + f20_month + 
                 f20_month_lag1 + f20_month_lag2 + arrest + 
                 release + as.factor(year) + as.factor(name_id), data = m)

stargazer(nb1, nb2, 
          covariate.labels = c("July 2015", "August 2015", "September 2015 (F20)", 
                               "October 2015", "November 2015", 
                               "Arrest", "Release"),
          title = "Negative Binomial Models of Online Searches and Media Coverage",
          omit = c("year", "name_id", "Constant"),
          label = "tab:nb",
          no.space = TRUE)

# Robustness for Dropping Countries ---------------------------------------

# Note: Core regression results in Table 3 of manuscript are robust to dropping any one country

# Drop countries
bc1 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
           any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[1]))
bc2 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[2]))
bc3 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[3]))
bc4 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[4]))
bc5 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[5]))
bc6 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[6]))
bc7 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[7]))
bc8 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[8]))
bc9 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[9]))
bc10 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[10]))
bc11 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[11]))
bc12 = glm(released ~ free20 + any_poc + any_torture + any_death_penalty +
            any_trial_legal + any_medical + as.factor(country), data = subset(c, country!=unique(c$country)[12]))

countries = unique(c$country)
lower = c(confint(bc1)[2,1], confint(bc2)[2,1], confint(bc3)[2,1], confint(bc4)[2,1], 
          confint(bc5)[2,1], confint(bc6)[2,1], confint(bc7)[2,1], confint(bc8)[2,1], 
          confint(bc9)[2,1], confint(bc10)[2,1], confint(bc11)[2,1], confint(bc12)[2,1])
upper = c(confint(bc1)[2,2], confint(bc2)[2,2], confint(bc3)[2,2], confint(bc4)[2,2], 
          confint(bc5)[2,2], confint(bc6)[2,2], confint(bc7)[2,2], confint(bc8)[2,2], 
          confint(bc9)[2,2], confint(bc10)[2,2], confint(bc11)[2,2], confint(bc12)[2,2])

as.data.frame(cbind(countries, round(lower, 3), round(upper,3)))
